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1. INTRODUCTION 


Consider the evolution partial differentiation equation u t = Lu, on a 
finite interval, where L is a hyperbolic operator. The solution u has a 
projection u on a finite subspace (which may for example consist of the 

first N modes in a Galerkin method, or N collocating points in the 

interval), and a numerical approximation u N generated by some spectral 

method. For linear operators it is known from the Lax equivalence theorem 
that if the scheme is consistent and stable, then u^ approximates P N u in 
some appropriate norm. If u is smooth, then the theorem implies that u N 
approximates the solution u in the same sense. 

In practice, one looks at the point values of u^ at the grid points and 
takes it as an approximation to the values of the true solution u at these 
points. We shall call this approach the realization of the computed solution 
via its grid-points value. The aims of the paper are: 1) demonstrate that 

when u is a complicated function, this realization will not produce 

acceptable results; 2) to suggest different ways for the realization of the 
solution in such cases. 

The following examples give a very clear illustration of the misleading 
results that may be obtained by pointwise realization. 

Example 1 

Consider the equation 




= u 


u(x,0) = u Q (x) 


0 < x < 2tt 


( 1 ) 



where u(x) and uq(x) are periodic functions and un(x) is a discontinuous 
function . If we expand uq(x) in Fourier series we get 


where 


u o (x) a E a k e 

k=-“ 


ikx 


(2a) 


a 


k 



(2b) 


The solution u(x) is thus given by 


u(x) = T' 


ikt ikx 
L* a k e e * 
k=-°° 


Suppose that (1) is solved numerically by the Fourier-Galerkin method, namely 
we seek a trigonometric polynomial of the form 


u N (x,t) 


£ b k (t)e 
k=-N 


ikx 


that satisfies 


9Un _ 9u^ 

3t 3x 



-N < k < N 


N 


u N (x,0) - •£ a k e 
k=-N 


ikx 


(3) 


From (3) it is clear that 


db k (t) 


ik b k (t), 


and 


dt 


-N < k < N 


(4) 



b k <0> ■ *k 


yielding the solution 

b k ( t ) - a k e . 

Therefore 

( \ ikt ikx / c \ 

u N (x,t) = 2-f a k e e t 5 ' 

Equation (5) implies that u N (x,t), obtained from the numerical solution (3), 
coincides with P N u(x,t), the Galerkin projection of u, thus yielding the 
best possible convergence of u N to P N u. However, since the Fourier series 
of u(x , t ) converges very slowly, the point values u N (xj»t) will not 
approximate well u(Xj,t). In general, one would witness the Gibbs phenomenon 
of overshoot in the neighborhood of the discontinuity and global oscillations 
all over the domain. In fact, even the initial approximation, u N (x,0), 
displays the same behavior in relation to uq(x). 

In the second example we show that the same phenomenon occurs even if the 
numerical initial point values do approximate the true initial point values to 
a high degree of accuracy. 


Example 2 

Consider the equation (1) where Uq(x) 


u Q (x,x) 


Ax 

A(2x - tt) 


is the saw-tooth function 


x < x 
x > x 


( 6 ) 


for some k, 0 < k < 2N-1, x = ~ (k+ V 2 )• • 



In the pseudospectral Fourier method we seek a trigonometric polynomial 


V N (x,t) 

V*’ 0 ■ s v t)eltx 

£=-N 


(7) 


such that 


3Vj i= 9 Vn 
3t 3x 


v N (x,0) = u Q (x,x). 


at the points x^ = , j = 0,...,2N-1 


(8a) 

(8b) 


Since Vjj is a polynomial of degree N, (8a) implies that 


9Vjj = 9Vn 
3t 3x 


(9) 


for all x (0,2 tt). Moreover, from (8b) it is clear that v^(x,0) is the 
(unique) trigonometric polynomial of order N that interpolates uq(x) at the 
points x j , j = 0,*.#,2N-1, thus 


N 

£ 

£=-N 


v N (x,°) = J) a A (x)e l£x = AF n (x,x) 


( 10 ) 


where 


a £ (x) 


2N-1 


-i£x. 


2N u 0 (x ;j> x)e J 

£ j=o 


(ID 


Performing (11) we get 


a Q (x) = A i [k - N + .5] 


( 12 ) 
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a t (x> 


= A 


1 - e 


-iir& 
N 


(k+1) 


2Nc 


1 - e 


-iiril 
N 


, , . it £ 

+ ict» 


- 1 


a * o 


(13) 


where 

C -N = °N = 2 ’ C £ = Ul * N * 

The numerical solution v^(x,t) of (9), (10) is 

v N (x,t) = v N (x + t,0) = AF n (x + t,x) (14) 

and upon manipulating (12), (13) one gets 

v N (x,t) = AF n (x,x - t) + At, (15) 

The trigonometric interpolant F^(x,x) collocates Uq(x,x) at the grid 
points x j • However, in between the grid points it oscillates. If we read 

the values of v^(x,t) at the grid points, then by (14) 

V N (x j’ t) = + t>X ) 

and unless t = . for some integer m, we will get solution that looks 

oscillatory. Thus, even though the initial approximation looks smooth at the 
grid points, when it evolves in time the oscillations will present themselves 
at the points Xj . 
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The conclusion one might draw from the above examples is that spectral 
methods (or any higher-order methods) are useless when applied to 
discontinuous function. A different approach is to look at a different 
realization of the numerical solution rather than the pointwise one. We will 
argue that high-order accurate information is contained in the numerical 
solution and demonstrate how that information can be extracted in such a way 
that accurate pointwise approximation to the true solution can be obtained. 


2. INFORMATION AND HOW TO EXTRACT IT 

Consider the linear equation 


u t = Lu 
u(0) = Uq 


(16) 


where L is a linear hyperbolic operator with variable coefficients and uq 
is a discontinuous function. For simplicity, we will restrict ourselves to a 
periodic, one (space) dimensional problem though the results are more general, 
(see Gottlieb and Tadmor [2]). Let v be the solution of the auxiliary 
problem 


v ■ - L v 
v(0) = v Q , 


(17) 


where vq is a C°° function. Because of the hyperbolicity of L, (17) is a 
well-posed problem. In Lemma 1 we quote the well-known Greenes identity. 



Lemma 1: Let u(t) 


and v(t) be the solutions of (16) and (17) at some 


level t, then 

(u(t), v(t) ) - (u 0 ,v Q ). (18) 

Assume now that (16) and (17) are discretized by the Fourier-Galerkin 
method. That is, we seek ujj and Vjj that are trigonometric polynomials of 
degree N such that for every k, |k| < N 


r 3u N T ikxN _ n 

' L v e ) ■ 0 

(19a) 

(u K (0) - u 0> e lkx ) - 0, 

(19b) 

(•“*■ ci? + v) - ° 

(19c) 

e ikx , (v N (0) - v Q ) = 0. 

(19d) 


We have also a Green identity for ujj and Vfj. 


Lemma 2: 

(u N (t),v N (t)) = (u N (0),v N (0)). (20) 

Proof : Since v N (t) and ujj(t) are N^-order trigonometric polynomials we 

use (19a) and (19c) to get 
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(; tT - L V v n ) " 0 

3v * 

(v Te“ + L v n> ■ °> 

and therefore 

It ( u s’ v n) ' (*"»•''») • ( u b> l * V ■ 0 

which implies (20). 

We will proceed by showing the relation of the RHS of (20) to that of 
(18). 

Lena 3: 

( U N (0) ’ V N (0) ) = ( u O ,V o) + e l (21) 

where 

lv n l 

Ui I < K — ^ (22) 

1 N S 

for every s. 

Proof ; From (19b) it is clear that 

(u N (0) - u Q , v N (0)) = 0. (23) 

Also, 

l( u 0» v N (0) ~ v q)I < K,u o" ,V N (0) ~ V 
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and since vq is a C function, 


" v n 11 s 

Jv (0) - v II < K — 

U N S 


for every s. 


(24) 


Now 

(V 0) ’ V N (0) ) = (u 0’V + ( U N (0) " V V N (0) ^ + ( U 0’ V N (0) " V 0 ) 
and in view of (23) and (24), 


(yo). V N (0) ) = (u 0’ v 0 ) + £ 1 

where 


and this proves the Lemma. 

From Lemmas 1 - 3 we can conclude: 

Theorem 1: Let u(t) and v(t) be the solutions of (16) and (17), 

respectively. Let ujj(t) and v N (t) be t ^ ie solutions of the Fourier- 
Galerkin approximations of (16) and (17). Then 

iv n i e 

|(u (t),v (t)) - (u(t),v(t))| < K — , for every s. (25) 

" N S 

The proof is an immediate consequence of (18), (20), and (21). 
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Assume now that the Fourier-Galerkin method described in (19c) and (19d) 
is stable, then vjj(t) approximates v(t) within spectral accuracy, that is 

II vB 

»v N (t) - v(t)» = e 2 < K— i . 

N 

We can, therefore, replace Vjj(t) In (25) and get 

(u N (t),v(t)) - (u(t),v(t)) + e 

where e Is spectrally small. We use now the fact that every C°° function 
v(t) can be obtained from some vq in (17). This is, in fact, one of the 
definitions of hyperbolicity. We can, therefore, state: 

Theorem 2: Let u(t) be the (nonsmooth) solution of (16) and let u N (t) be 

the solution of the spectral Galerkin approximation to (16), Then for any C°° 
function v(t) 


(u N (t),v(t)) - (u(t),v(t)) + e (26) 
where e is spectrally small. 

Thus, u^(t) approximates weakly u(t) within spectral accuracy. It is 
in this sense that ujj(t) contains a highly accurate information about 
u(t). We will show later how to use this information in order to obtain 
spectral accurate approximation to the grid-point values of u(t). 
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We turn now to the pseudospectral Fourier case. Here we need some 
preprocessing of the initial data in order to prove the same result as in 
Theorem 2. 


Theorem 3: 

satisfies 


Let u^(x,t) be a trigonometric polynomial of order 



Lu 


N 


at x 



0 , 


2N-1 


N 


(^(0) - u 0 , e ikx ) = 0, 


k| < N, 


that 


(27) 


(i.e., Ujj(x,t) is the solution of the pseudospectral Fourier scheme, but 
initially ujj(x,0) is obtained by the Galerkin projection). 

Then for every smooth function u(x,t) 


2N-1 2tt 

F £ u«(x ,t) v(x ,t) = J u(x,t) v(x,t)dx + e 
w j =0 3 3 0 


(28) 


where e is spectrally small, provided that the pseudospectral approximation 
is stable. 


Proof : Let v ^ be the solution of the pseudospectral Fourier approximation 

of (17a) and let v^(0) be the Galerkin projection of vq, that is 


( V N (0 > " V e±kX ) = °* l k l < N * 


( 29 ) 


From (27) and the analog equation for v N> one gets 
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2N-1 

S’ u N (lt J ,c) Vv 0 


JL 

N 


N 


£ 

j=0 


u N (x j ,0)v N (x j’ 0) * 


(30) 


From the exactness of the trapezoidal rule for polynomials of degree 2N, we 
conclude 



Note that the initial functions v N (x,0) and u N (x,0) are not the 
interpolants of uq and vg as in the usual pseudospectral methods but rather 
the Galerkin approximation to ug and vg. We recall now Lemma 3 and 
equality (18) to establish (28). The proof is thus completed. 

It is interesting to note the way in which the information is contained. 
The interpolant of ug looks smooth at the grid points, whereas the Galerkin 
approximation of ug looks oscillatory on the grid points. It means that in 
order to preserve the information one has to require initially oscillatory- 
looking solution. The information is preserved m the structure of the 
oscillations. 

We will show now a way of using (26) and (28) in order to construct a 
better approximation to u(xj,t) then the one given by u^(xj>t) (here 
u N (x,t) is given by either the Galerkin method or the pseudospectral method). 

From (28) and (26) it is clear that in order to get a good approximation 
to u(y,t) at some point y (0,2tt), we need to find a function v^(x,t) 


2tt 


f u(x,t) 


Vy(x, t)dx = u(y,t) + Ej, 


such that 
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where is spectrally small. By (26) we will have 


2ir 

J u^(x,t) v (x,t)dx = u(y,t) 
n y 


+ e + e. 


for the Galerkin method and 


1L 

N 


2N-1 

E 

j=0 


^(x^t) v y (x^,t) = u(y,t) + e + e 


for the pseudospectral method. 

For conveniency we will shift the interval [0,2ir] to [— 
p(x) be a C°°-function vanishing outside the interval 


p(0) = 1. 


Let D p (x) be the Dirichlet kernel, namely 


D p (y> 


l 

2* 


E^ 

|k|<p 


1 sin(p+ V 2 ) y 

TH TuTTyTT) 


We set now 

v y (x) = t|> 0,P (x) = 0 1 p(0 ^>0^(0 1 y). 


One can prove (see [2 ]) that 


J u(x)i|> 0 ,P (y— x)dx = u(y) + 


where is spectrally small. 

Thus, it is possible to extract accurate pointwise values 


(32) 

(33) 

,ir]. Let 
satisfying 

(34) 

(35) 

(36) 

(37) 

from u N (x). 
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3. NUMERICAL RESULTS 

In this section we demonstrate the efficacy of the smoothing procedure 
outlined above. As a test function we have chosen the piecewise C -f unction 


f(x) = 


X 

sln *2 


-sin j 


0 < X < 7T 


TT < X < 2tr. 


Denote its spectral approximation by f^(x), and let ?^(x) be the 
pseudospectral approximation to f(x). It is evident from the first column of 

A 

Tables I and III that ~ the spectral approximation sampled at 

y v = vtt/N - do not approximate f(y^) within spectral accuracy. In fact, 

A A 

the error committed by f 100 (y ) is only half of that committed by f,,(y )• 

12o v '64 

Regarding the pseudospectral approximation, f^(x), it, of course, collocates 

the exact values at the sampling grid points, ?^(y^) = f(y^); yet, in between 

these gridpoints, ^ N (y v+ = + V 2 H/N) approximate f(y v+ l^) within 

first-order accuracy only, as shown in the first column of Tables II and IV. 

In order to construct our regularization kernel in (36), we define the 

cut-off function p(g) = p (?) to be 

a 


P a ( 5) = 


exp 

r - 1 


Is I < 1 


otherwise 


namely, p^(£) is a C - function whose support is the interval l^j < 1. 
ip to be used in (36) is of the form 


* 0 »P(y) =_L 0 (0 -1 y) Sill (p+ ]b )y/6 

* vy; 2x0 < V° y) sin y/20 


( 40 ) 
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The post-processing procedure of the spectral approximation 
* 0 D 

convoluting against namely 


involves 


f (x) 


2i? f ?,WP(V) 


sin (p+ V2 )(x-y)/0 
sin (x-y)/20 dy 


(41) 


where x is a fixed point of interest. (In practice, we use the trapezoidal 
rule to evaluate the right-hand-side of (41) taking a large number of 
quadrature points.) 

The parameter 0 was chosen as 


0 = IT* | X - TT | ; 


(42) 


this guarantees that ij; is so localized that it does not interact with 
regions of discontinuity. 

It should be noted, in this stage, that if 0 was so chosen to be the 
same for each x, (and not as in (42)), the formula (41) admits a simpler 
form; that is, if 

i|» 0,p (y) = S a v eiky (43) 

• k=-» 

then 

f(x) ~ V f(k)a, e . (44) 

k=-N 

This procedure can be carried out efficiently in the Fourier space. 

Next, we turn to the post-processing for the pseudospectral approximation 

A 

f^(x) which is simpler than (41). In fact, in this case 



-16- 


9 2N " 1 n « 

f U) ~ ^ £ f(y v H 0 ,p ( x ~y v )* (45) 

v=0 

Note that carrying out the smoothing procedure defined in (45) does not 
involve any extra evaluation of f(y) in points other than y^, in contrast 
to spectral smoothing procedure in (41). As before, the parameter 6 was 
chosen according to (42). We have yet to determine the parameters p and 
a. The parameter p must be equal to N^ for 0 < $ < 1, in order to 

assure infinite accuracy, (In our computations, 3 « .8.) Finally, we feel 
that a is problem dependent and we chose a = 10. We have not tuned the 
parameters to get optimal results; further tuning may improve the quality of 
our filtering procedure. 

In Tables I, II, III, and IV we give the results of the smoothing 
procedure at several points in the domain. The pointwise values are now 

recovered with high accuracy. The first column in each table indicates the 
points in which the procedure was performed. We limited ourselves to four 
points in the interval (0,ir) because of the symmetry of the function f(x). 

A 

The second column gives either the spectral approximation f^(x) or tlie 
pseudospectral approximation ?^(x), N = 128 in Table I and II and N - 64 in 
Tables III and IV. The third column gives the smoothed results, when filtered 
by (41) on (45), at the same points as in column I. 

The results indicate the dramatic improvement obtained by the smoothing 
procedure. Moreover, note that the error committed by ? 100 (or f Q ) is 

J.ZO I/O 

A 

better than the one committed by f^ (or f^) only by a factor of 2 whereas 
after the post-processing the error improves by a factor of 10 4 . 
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Table I. Results of smoothing of the spectral 
approximation of f(x), N = 128 


TTV 

X = o" 

v 8 

v equals 

|f(x v ) - f N (x v )| 

1 * * l 

l £ - £ n *1 

at x - * v 

2 

3.2 (-3) 

5.8 (-10) 

3 

5.2 (-3) 

7.9 (-10) 

4 

7.8 (-3) 

6.3 (-10) 

5 

1.1 (-2) 

1.1 (-10) 


Table II. 

Same as Table I for the pseudo- 
spectral approximation ? N (x)* 


*v+l/2 ■| <v+l/ 2 ) 

v equals 

|f(% + v 2 ) ' Vvv 2 }| 

l £ - *. * *1 

at X ‘ V 1/2 

2 

5 (-3) 

7 (-10) 

3 

8.1 (-3) 

7.9 (-10) 

4 

1.2 (-2) 

6.4 (-10) 

5 

1.8 (-2) 

1.2 (-10) 
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Table III. Results of smoothing of the spectral 
approximation of f(x), N = 64 


TTV 

x - — 5- 
v 8 

v equals 

|f(x v ) - f N (x v )| 

1 A * ■ 
l f ‘ f N *1 

at x = x 

V 

2 

6.4 (-3) 

4.8 (-6) 

3 

1 (-2) 

5.9 (-6) 

4 

1.5 (-2) 

7.7 (-6) 

5 

2.3 (-2) 

8.9 (-6) 


Table IV. 

Same as Table III for the pseudo- 
spectral approximation, ?^(x). 



V l/ 2 "if <V+1/2) 

l f(x v+V 2 ) " *, ( WI 

If 

- *■ * *1 

v equals 


at 

x = V 1/2 

2 

1 (-2) 


4.1 (-6) 

3 

1.6 (-2) 


6 (-6) 

4 

2.4 (-2) 


7.8 (-6) 

5 

3.6 (-2) 


8.9 (-6) 
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4. A DIFFERENT METHOD FOE EXTRACTING INFORMATION 

In this section we would like to present a different approach for 
extracting the information from an oscillatory solution. The idea is to 
subtract from the solution those oscillations that correspond to the saw-tooth 
function discussed in Example 2. This leads to the following procedure: 

» * ilx 

Let u„(x,t) = u e , be the solution of the pseudospectral 

N l - -N 1 

approximation to a hyperbolic problem. We try to find an unknown smooth 
function and a (oscillatory) saw-tooth function F N (x-t,x s ) with an unknown 
jump 2trA at an unknown location x g such that 


2N-1 


H 


2 U N (x j’° ' AF N (x j* X s ) “ C " 2 \ 6 


ilk 


1 2 


(46) 


j“0 


P 

1*0 


is minimized. Note that we have 2p + 3 unknowns in (46): A, x s » c and 

2P values of b (i * 0) . 

X 


The conditions for local minima of H are found from the following 
2p + 3 equations: 


3H 

3A 


0 


2N-1 



u . F - AFT - cF. 
j j 1 J 


F. 

J 




= 0 


(47) 


where F^ = F N (xj,x g ), u^ 


u^ (Xj » t ) . Also, 
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3H 

3c 


2N-1 

0 "> £ 

J=0 


Ux. 


u. - AF - c - V b, e J , = 0 
J J 2^ * 


£=-p 
I 1*0 


2N-1 


IS. = o ==> V f: u. - af: f. - cf: - f: V 

3s LmJ J J J J J J Z-f 


b rt e 


iX j A 


j=0 


*=-p 

1*0 


where FT 
J 


= 3F N ( x., x s )/ 3 s = ^ 


3a. (s) i£x 


3s 


e ^ ; and 


A— N 


___ = o ==> b = u - Aa, m = 1,2, ... ,p 
3b m m m’ 1 1 ’ ’ 

m 


where 


2N-1 


u = 


_ L V 

2Nc 

m . ^ 


-i£x. 


u N (Xj)e 


j=0 


Substituting (50) into (47), (48), and (49) we get, respectively: 


- Aa 0 - 


c = 0 


£ 

u I>p 





1 . u. 
-I A 



-A l 


= 0 


) 


where a'(s) = 3a^(x)/3s. 


Next, we combine (52) and (53) to get 


(48) 

(49) 

(50) 


(51) 

(52) 

(53) 
single 
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nonlinear equation for s: 


2X a - f . E c 


H c o a — 0 U fl 2 C 


(54) 


where all sums run over p < |&| < N. 

Equation (54) is solved iteratively for s. Having found s, one 

immediately obtains from Example 2 all the a (s)'s. Then from (50) we have 

jc 

the b m 's, and A from (52). Finally, having A we find c from (51). 

The minimum thus obtained may be a local one while we are seeking a global 
minimum. This means that in practice one searches for the global minimum. 

We now give an example that illustrates the efficacy of the procedure. We 
solve the following problem: 


^ U N ** U N 

— _ + = o 

3t 3x ’ 


0 < x < 2 tt, t > 0 


(55) 


u N (x,0) = 


sin y 


4 X 
" sin 2 


0 < X < 71 


TT < X < 2tT 


(56) 


u N (0,t) = U N (2TT,t). 


(57) 


We ran the problem on several grids and exhibit here the numerical results for 
the case N = 8 (i.e., 16 subintervals in the domain (0,2ir)). The 
unadulterated results at t = ir/2N on the grid points are shown in Figure 1. 



exact solution 


x unsmoothed pseudo- 

spectral solution 

o smoothed solution 

(N = 8) 



X 


Figure 1 



Table V 


j 

exact 

solution 

error 1 = 

| exact-unsmoothed | 

error 2 = 

| exact-smoothed | 

error 1 
error 2 

0 

9.80 

X 

10- 2 

5.86 

X 

10“ 5 

5.86 

X 

io - 5 

1.00 

1 

9.80 

X 

10 -2 

1.24 

X 

10 -2 

5.86 

X 

10~ 5 

211 

2 

2.90 

X 

10 -1 

2.57 

X 

10 -2 

6.30 

X 

10 -5 

408 

3 

4.71 

X 

10" 1 

4.13 

X 

10“ 2 

7.33 

X 

IO -5 

563 

4 

6.34 

X 

10' 1 

6.15 

X 

10“ 2 

9.30 

X 

10 -5 

661 

5 

7.73 

X 

10 _1 

9.11 

X 

IO -2 

1.31 

X 

10 -4 

695 

6 

8.82 

X 

10 -1 

1.43 

X 

10 -1 

2.16 

X 

10 -4 

662 

7 

9.57 

X 

10 -1 

2.70 

X 

10" 1 

4.42 

X 

10' 4 

611 

8 

-9.95 

X 

10 -1 

1.00 

X 

10° 

1.10 

X 

10" 2 

91 

9 

-9.95 

X 

10 -1 

2.68 

X 

io -1 

1.34 

X 

io -3 

200 

10 

-9.57 

X 

10 -1 

1.42 

X 

10" 1 

4.42 

X 

IO -4 

321 

11 

-8.82 

X 

10 -1 

9.07 

X 

io -2 

2.16 

X 

10 -4 

420 

12 

-7.73 

X 

10 -1 

6.12 

X 

io -2 

1.32 

X 

10 -4 

464 

13 

-6.34 

X 

10 _1 

4.11 

X 

IO -2 

9.30 

X 

10" 5 

442 

14 

-4.71 

X 

10 _1 

2.55 

X 

10 -2 

7.32 

X 

io - 5 

348 

15 

-2.90 

X 

10 -1 

1.22 

X 

10 -2 

6.30 

X 

io - 5 

194 


We then post-processed these u N (x^ ,tt/ 2N) values according to the procedure 
described above. The filtered values are shown on the same graph, and the 
errors listed in Table V are computed before and after processing. The 
dramatic improvement is evident. 

Next we demonstrate the procedure in the case of the Euler equation for 


gas dynamics. Because the physical problem involves inflow, outflow, and no- 
flow boundary conditions, periodicity could not be imposed and we use the 
Tchebyshev, rather than Fourier, pseudospectral method. 
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The physical problem is that of a wedge, inserted as a zero angle of 
attack, into a uniform supersonic flow of an ideal gas with y = 1.4. An 
oblique shock develops in time and the flow reaches, after a while, a steady 
state. The time— dependent Euler equations in two-space dimensions were 
discretized by the pseudospectral Tchebyshev method in space with an 8x8 grid 
and a modified Euler scheme was used for the time discretization. Since we 
are interested in the steady state only, the accuracy for the time integration 
is of little importance. In order to be sure that a steady state is reached, 
the code was run until all physical quantities did not change to 11 
significant figures over a span of 100 time steps. The values of the density 
in the steady state at the grid points together with the grid points 
themselves are given in Table VI. 


Table VI. 


p 

Y 


1 .862 

1.851 

1.869 

1.871 

1.837 

1.865 

1 .892 

1.885 

1.878 

1 . 


1 .862 

1.870 

1.867 

1.820 

1.870 

1.954 

1.899 

1.803 

1.759 

.961 


1.862 

1.854 

1.852 

1.904 

1.877 

1.770 

1.782 

1.864 

1.900 

.853 


1.862 

1.871 

1.876 

1.812 

1.838 

1.969 

1.975 

1.884 

1.841 

.691 


1.862 

1.848 

1.842 

1.935 

1.899 

1.703 

1.710 

1 .890 

1.984 

.5 


1.862 

1.883 

1.894 

1.729 

1.832 

2.429 

2.994 

3.255 

3.316 

.308 


1.862 

1.808 

1.810 

2.387 

3.133 

3.375 

3.224 

3.054 

3.002 

.146 


1.862 

2.115 

2.868 

3.288 

3.176 

2.965 

3.006 

3.136 

3.187 

.038 


1.862 

3.083 

3.046 

2.975 

3.087 

3.108 

3.024 

3.013 

3.016 

0 

X 

0 

.038 

.146 

.308 

.5 

.691 

.853 

.961 

1 . 
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Note that the raw data in Table VI seems to indicate roughly the same y-shock 
location at Xq = 1, xj = .961, and X 2 = .853, namely between the grid points 
y^ = .3086 and y$ - .500. This means that because of the coarse Tchebyshev 
grid the shock location cannot be resolved to better than 20% of the domain. 
In fact, the correct shock locations at those x-stations are y = .434 for 
x 0 , y = *417 for x L and y = .370 for X£. 

In the present case it is not necessry to employ a saw-tooth piecewise 
smooth function, as was done in the previous section, because there is no need 
to preserve periodicity. Instead, we subtract from the oscillatory data an 
expansion of the Heaviside function, S(y,y g ): 


S( y,y s ) 


j d l + d 2 

k 


-i < y < y s 
y s < y < 1 


(58) 


where d^ 9 the state ahead of the shock, and d 2 , the magnitude of the 
discontinuity, are constant. The description here of S(y,y s )> as if 
independent of x, has to do with the fact that the two-dimensional results of 
the pseudospectral algorithm were post-processed at fixed x-stations. The 
expansion of S(y,y g ) is given by 

s N (y,y s ) = 2 A £ (s)T £ (y) (59) 


where T^(y) is the Tchebyshev polynomial of order 
T^(y) = cos[fc cos -1 (y)], and 
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A Q (s) - (s + j)/N 

A^s) * sin[^|j- (s + j )] /N sin ; 1 < i < N-l 

A^s) - sin[(s + j)]/2N. 


If s 


is an integer, 


then on the grid points. 


y j 


cos(irj/N) . 


S N (y j» y s ) = S(y j» y s ) * 


(60) 


The which we wish to minimize is now, at any given x-station: 






p<N 

- E 

t = i 


b „ T » 

i i 


<vi : 


(61) 


( 1 1 < j < N-l 

(2 j = 0, N 


(62) 


Differentiating (61) with respect to the parameters dj, d 2 > s and 
b(l<A<p<N), using the orthogonality relations for the Tchebyshev 
polynomials and manipulations similar to those used in the previous section, 
we get p + 3 nonlinear algebraic equations which are completely analogous to 
(50) - (53). They are: 


, £-l,2,...,p. (63) 

p 0 - d 2 ^ - d l = 0 


(64) 
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N 

E 


N 


C A \ P £ " d 2 E C £ \ = ° 

£=p+i * * * z * 11 


(65) 


where 


N 

£=^j+l 


N 


E c « a « p 0 ■ d o y c. a, a: = o 

* i I t 2 .iU 


£=p+l 


III 


N 


»t ’ TC7 ScT^WV 


4 j=0 j 


A' -fcV-). 


(66) 


(67) 


( 68 ) 


Again, we combine (65) and (66) into a single nonlinear equation for the shock 
location index, s: 


y. A' p Vc - 5^0. A p yc A A' = 0 (69) 

^ Z Z y Z ^ Z Z ^ZZ^Z^ZZZ 

where all the sums are from Z = p+1 to £ = N. 

The procedure for extracting the shock location, jump magnitude and smooth 
part of the solution from the raw data p(x,y^) (given in Table VI) is 

exactly the same as described above for the Fourier problem. 

For the wedge-flow problem considered here, this procedure applied in the 
case of a coarse net (N = 8), located the shock with an error only in the 
fourth significant figure. The smooth part was recovered to within 1% at the 
worst field point. 
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Conclusion 

We have demonstrated that the realization of a numerical solution via its 
grid-point value may be misleading when the true solution has a complicated 
structure which is not resolved by the grid. We have shown, however, that the 
numerical solution does contain highly accurate information about the solution 
and we suggested two ways of extracting this information. 

The analysis outlined in this chapter is a linear one (though the 
procedure was applied also to nonlinear problems). However, in [28] Lax has 
argued that more information about the solution is contained in high 
resolution schemes even in the nonlinear case. In fact, using notions from 
information theory. Lax has shown that the e-capacity of the set of 
approximate solutions is closer to the e-capacity of the set of the 
projections of exact solutions if the numerical scheme is a high-order scheme. 

In the area of digital filters one always processes the data in order to 
overcome the Gibbs phenomenon. If we look at the initial conditions as an 
input signal and at the final result as the output signal, the idea of 
filtering is a natural one. 
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